/* GNU Mailutils -- a suite of utilities for electronic mail
   Copyright (C) 1999, 2000, 2001, 2005, 2007 Free Software Foundation, Inc.

   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 3 of the License, or (at your option) any later version.

   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General
   Public License along with this library; if not, write to the
   Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
   Boston, MA 02110-1301 USA */


/**
 * \file silhouette.cc
 * \author Ice and Magic
 * \date 20 Avril 2011
 * \brief Machine abstraite pour le concept de "Silhouette"
 */

#include <iostream>
#include <sstream>
#include "silhouette.h"
#include "global.h"


using namespace std;
/**
 * \brief matrice contenant les puissances $X_i^k$
 * $X_i$ un chiffre d'un entier $X$
 * $k$ un entier
 *
 * Le calcul des puissances n'est fait qu'une seule fois. On peut
 * faire mieux, en les codant en dure ... ou alors mieux, de les
 * calculer sur ma machine cible puis de les inclure avant de 
 * compiler le programme qui fait l'exploration des silhouettes.
 * TODO: je vais peut être le faire avant le release final
 */

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////

Silhouette::Silhouette(__type t){
#pragma ivdep
for(int i=0;i<__DIGITS;i++)
	digit_occurrences[i]=0;	

for(;t>0;t/=10)
	(*this)[t%10]++;
}
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
/*
inline Silhouette&
Silhouette::operator=(Silhouette& s){
#pragma ivdep
for(int i=0;i<__DIGITS;i++)
	digit_occurrences[i]=s.digit_occurrences[i];
return *this;

}
*/
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////

inline int
Silhouette::nb_digits(){
int res=0;
#pragma unroll(10)
for(int i=0;i<__DIGITS;i++)
	res+=digit_occurrences[i];
return res;
}
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
Silhouette::~Silhouette(){
/*Don't bother me again, i'm doing nothing*/
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////

inline bool
Silhouette::operator==(Silhouette&	s){

// le coef du Zero n'a aucun impact
bool res=true;

#pragma unroll(9)
for(int i=1;i<__DIGITS;i++)			
	res &= (s.digit_occurrences [i] == digit_occurrences [i]) ;

return res;
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
inline short &
Silhouette::operator[](int i) {
return	digit_occurrences[i];
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
std::string	
Silhouette::humain_friendly_print(){
stringstream res;
res<<"<";
for(int i=__DIGITS-1;i>0;i--)
	for(int j=0;j<(*this)[i];j++)
		res<<i;
res<<">";
return res.str();
}
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
ostream&	
operator<<(ostream& out,  Silhouette& s){
out<<"{";
for(int i=0; i<__DIGITS; i++){
	out<<s[i];
	if(i<__DIGITS-1)
		out<<", ";
	}
out<<"}";
out<<s.humain_friendly_print();
return out;
}
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////

__type	
Silhouette::compute_F_x(int k){
__type result;
		// TODO: virer cette fonction de la version release.
// coef* 1 ^ k = coef
// coef* 0 ^ k = coef
result=(*this)[1];

for(int i=2;i<__DIGITS;i++)	
	// $\sum{i}{(coef_{X_i}*{X_i}^k)}$
	result+=(*this)[i]*MyPower(i,k);	

return result;
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////

__type	
Silhouette::compute_F_x(const __type* powers){
__type result;

// coef* 1 ^ k = coef
result=(*this)[1];	


//__type	ir[__DIGITS];

#pragma unroll(8)
for(int i=0;i<__DIGITS-2;i++)
	result+=this->digit_occurrences[i+2]*powers[i+2];

//#pragma unroll(8)
//for(int i=2;i<__DIGITS;i++)
	// somme_i(coef_Xi* Xi^k)
//	result+=ir[i];				//(*this)[i]*powers[i];	

return result;
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////

bool
Silhouette::is_digital_invariant(const __type* power_k){
Silhouette	f_x(compute_F_x(power_k));
return (*this)==f_x ;
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////

/**
 *  complexité: k (nombre maximal de chiffres) * 9 (nombre de chiffres)
 *  O(k)
 *
 */
set<__type>
Silhouette::collect_invariants(int k){
set<__type> result;

#ifdef	FULL_VERBOSE
cout<<endl<<"\tDEBUG : Silhouette::collect_invariants ..... k = "<<k<<"\n\t Liste des F(x) Générés";
#endif



//TODO prefetcher power_matrix

// cette silhouette peut être celle un point fixe d'ordre k+3 pour une __CAPACITY=19
// voir preuve \cite ref3 ref1


//#pargma prefetch(power_matrix)
for(int i=__MAX(k,1); i<=__MIN(k+1,__CAPACITY-1);i++){
	
	__type*	powers;

	powers=power_matrix[i];


	//TODO prefetcher powers
	// (*this) = X, calculer F(X)
	__type	f_x;
	#pragma prefetch(powers)
	f_x=compute_F_x(powers);

	// Silhouette de F(x)
	Silhouette F_x(f_x);

	// Calculer F(F(X))
	__type	f2_x=F_x.compute_F_x(powers);
	Silhouette F2_x(f2_x);

	#ifdef	FULL_VERBOSE
	cout<<"\n\t PUISSANCE i = "<<i<< "\tF (x) = "<< F_x<<"\t NB = "<<f_x<<" f2_x = "<<f2_x;
	#endif

	// F(F(X)) = F(X)	==> Point fixe d'ordre k potentiel
	if(F_x == F2_x ){
		// F(F(X)) a k chiffres?
		if(F2_x.nb_digits() == i){
			// Oui, X est un attracteur et  F(F(X)) est un point fixe.
			result.insert(f2_x);
			#ifdef	FULL_VERBOSE
			cout<<"\t Point fixe =  "<<f2_x;
			#endif
			}
		}
	else	{

		/* Cycle potentiel => éviter d'y entrer
		 *Si cette Silhouette est attractrice, il y aura une autre
		 * Silhouette qui convergera vers le même point fixe.
	 	 */
		}
	}
return result;
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////


/**
 * Cette fonction fait, essentiellement, l'exploration des éléments de la classe d'équivalence
 * des entiers définie par notre relation. L'exploration se fait par génération récursive de ces
 * éléments. Ce qui crée un arbre implicite d'ordre 9 (9 chiffres significatifs 1..9).
 *
 * \param nb un élément de la classe d'équivalence. Util
 * \param s sa silhouette.
 * \param remained_digits profondeur du sous arbre.
 * \param nb_digits profondeur de tout l'arbre. Au départ nb_digits = remained_digits.
 * \param upper_bound 
 * \return l'ensemble des points fixes collectés.
 *
 * Elle fait aussi, étant donné que l'exploration est fait, le test de point invariant pour chaque
 * noeud de l'arbre. Il ne reste qu'à les sélectionner.
 */


//TODO: le paramètre nb est a virer de la version release

set<__type>
generateur_silhouettes(
			#ifdef FULL_VERBOSE
			__type nb,
			#endif
			Silhouette s, int upper_bound, int remained_digits,int nb_digits){
						// Il est révolu le temps où on disait que la
						// la récursivité coûtait chère.
						// profondeur maximale = 18
						//
						// Edit Bon, étant les perfs publiées sur le forum,
						// je vais peut retirer cela ;)
						// J'aurais dû penser à une solution itérative ;)


// Une feuille de l'arbre d
if(remained_digits == 0){
	#ifdef	FULL_VERBOSE
	cout<<endl<<"DEBUG : generateur_silhouettes"<< endl<< " S = "<<s << "\t NB = "<<nb;	
	#endif
	return s.collect_invariants(nb_digits);
	set<__type> res;
	return res;
	}

// Ce n'est pas une feuille.

set<__type> 	result;
set<__type>	results[upper_bound];


	for(int i=1;i<upper_bound;i++){
		Silhouette si=s;
		si[i]++;
		results[i]=  generateur_silhouettes(
				#ifdef FULL_VERBOSE
				nb*10+i,
				#endif
				si,i+1,remained_digits-1,nb_digits);	
		}
	
	

// Collect des tous les points fixes des fils
for(int i=1;i<upper_bound;i++)
	result.insert(results[i].begin(),results[i].end());

// le noeud courant, est aussi un élément de la classe d'équivalence!
// ==> calculer les ses points fixes s'il est attracteur.
set<__type>	sp_set=s.collect_invariants(nb_digits-remained_digits);	

// Collection finale
result.insert(sp_set.begin(), sp_set.end());

return result;
}

//////////////////////////////////////////////////////
//////////////////////////////////////////////////////


